#include "ICritAlgo.h"

#include "TrajSimul.h"
#include "SwitchingCheck.h"

#include <cmath>
#include <iostream>

ICritAlgo::ICritAlgo()
{

}

ICritAlgo::~ICritAlgo()
{

}

double ICritAlgo::solve(Clock clock,Output output,Pulse *field,Pulse *current,MultiLayer multiLayer,NewParameters params)
{

  int run=params.getNrun();
  if(params.getTemperature()==0)
    run=1;
  double proba=0,nrun=0,nswitch=0;
  double iCritmax=params.getCurrent()["max"];
  double iCritmin=params.getCurrent()["min"];
  double iStep=params.getCurrent()["step"];
  MultiLayer result=multiLayer;
  bool bmin,bmax;
  current->setLevel(iCritmin);
  for(nrun=0;nrun<run;nrun++)
    { 
      result=TrajSimul().solve(clock,output,field,current,multiLayer,params);
      bmin=SwitchingCheck(result[multiLayer.getSolveLayer()],params.getTorque()["tolerance"]).getCheck();
      if(bmin) 
	{
	  nswitch++;
	  std::cout<<"switching"<<std::endl;
	}
      else std::cout<<"not switching"<<std::endl;
    }
  proba=nswitch/run;
  if(proba>=0.5)
    bmin=true;
  else bmin=false;
  //std::cout<<"switching"<<std::endl;
  //else std::cout<<"not switching"<<std::endl;
  //std::cout<<"imin :"<<proba<<std::endl;
  nswitch=0;
  current->setLevel(iCritmax);
  for(nrun=0;nrun<run;nrun++)
    {
      result=TrajSimul().solve(clock,output,field,current,multiLayer,params);
      bmax=SwitchingCheck(result[multiLayer.getSolveLayer()],params.getTorque()["tolerance"]).getCheck();
      if(bmax) 
	{
	  nswitch++;
	  std::cout<<"switching"<<std::endl;
	}
      else std::cout<<"not switching"<<std::endl;
    }
  proba=nswitch/run;
  if(proba>=0.5)
    bmax=true;
  else bmax=false;
  //std::cout<<"imin :"<<proba<<std::endl; 
 if(!bmin||bmax)
    {
      std::cout<<"bmax ou !bmin"<<std::endl;
      return 0;
    }
 while(fabs(iCritmax-iCritmin)>iStep)
    {
      bool btemp;
      double iTemp=(iCritmax+iCritmin)/2;
      current->setLevel(iTemp);
      nswitch=0;
      for(nrun=0;nrun<run;nrun++)
	{
	  result=TrajSimul().solve(clock,output,field,current,multiLayer,params);
	  btemp=SwitchingCheck(result[multiLayer.getSolveLayer()],params.getTorque()["tolerance"]).getCheck();
	  if(btemp) nswitch++;
	}
      proba=nswitch/run;
      if(proba>=0.5)
	btemp=true;
      else btemp=false;
      if(btemp&&bmin)
	{	
	  iCritmin=iTemp;
	  bmin=true;
	}
      if(!btemp&&!bmax)
	{	
	  iCritmax=iTemp;
	  bmax=false;
	}
    }
  return iCritmin;
}
